#include "fox_mul.h"
#include <unistd.h>
#include <sys/time.h>

//Multiplication sequentielle
void block_mult(int * A,int* B,int* C,int taille_block)
{
  int i,j,k,n;
  int N = taille_block*taille_block;
  for (n = 0 ; n < N ; n++){
    for(k = 0 ; k < taille_block ; k++){
      i=n/taille_block;
      j=n%taille_block;
      C[n] += A[i*taille_block+k]*B[k*taille_block+j];
    }
  }
}

//taille_mat : nbre d'element sur une ligne.
int * parseur(FILE*fd,int * taille_mat, int * nb_block)
{
  //On récupere le nombre de block.
  char bff[3]; //assuming nbr of blocks < 99
  fgets(bff,4,fd);
  int l = strlen(bff);
  bff[l-1] = '\0';
  *nb_block = atoi(bff);
  
  
  //Calcul de la taille de la matrice.
  
  int nb_ligne = 0;
  int c;
  do {
    c = fgetc (fd);
    if (c == '\n') 
      nb_ligne++;
  } while (c != EOF);
  //printf ("The file contains %d lines\n",nb_ligne);
  rewind(fd);
  
  *taille_mat = nb_ligne;
  
  int * matrice = malloc(sizeof(int)*(*taille_mat)*(*taille_mat));
  int k = 0;
  char buf[100];
  fscanf(fd,"%s",buf);
  while(-1 != fscanf(fd,"%s",buf))
    {
      matrice[k] = atoi(buf);
      ++k;
    }
/*   ////////////////// */
/*   int i=0,j; */
/*   while(i<taille_mat*taille_mat) */
/*     { */
/*       for(j=0 ; j<taille_mat ; ++j) */
/* 	{ */
/* 	  printf("%d ",matrice[i]); */
/* 	  i++; */
/* 	} */
/*       printf("\n"); */
/*     } */

  //////////////// 
  return matrice;
}

void print_block(int * block, int taille_block)
{
  int i,j;
  for (i=0;i<taille_block;i++)
    {
      for (j=0;j<taille_block;j++)
	{
	  printf("%d ",block[i*taille_block+j]);
	}
      printf("\n");
    }
  printf("\n");
}

int main(int argc, char* argv[])
{
  if (argc < 3)
    {
      printf("utilisation : %s matA, matB, res.data",argv[0]); 
    }
  int tag = 42;
  int myrank;
  int size;
  int nb_block;
  int taille_mat;
  int taille_block;
  int block_s;
  int buf[2];
  MPI_Init(NULL,NULL);
  MPI_Status status;
  MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
  MPI_Comm GridComm;
  MPI_Comm_size(MPI_COMM_WORLD,&size);
  int * A;
  int * B;
  struct timeval init,mid,end; 
  if(myrank == 0)
    {
      gettimeofday(&init,NULL);
      char * matA = argv[1];
      char * matB = argv[2];
      //Chargement des matrices.
      //Matrice A:
      int taille_matA;
      int nb_blockA;
      FILE * fdA = fopen(matA,"r");
      A = parseur(fdA,&taille_matA,&nb_blockA);
      fclose(fdA);
     
      //Matrice B:
      int taille_matB;
      int nb_blockB;
      FILE * fdB = fopen(matB,"r");
      B = parseur(fdB,&taille_matB,&nb_blockB);
      fclose(fdB);
      //Verification de la cohérence des matrices
      if (taille_matB != taille_matA){
	printf("Les matrices sont de tailles différentes.\n");
	return 0;
      }
      if (nb_blockB != nb_blockA){
	printf("Les blocks sont de tailles différentes.\n");
	return 0;
      } 
      if(taille_matA%nb_blockA != 0)
	{
	  printf("Le nombre de blocks(%d) ne divise pas la taille de la matrice(%d)\n",nb_block,taille_mat);
	  return 0;    
	}
      if(size != nb_blockA*nb_blockA){
	printf("Processus trop ou pas assez nombreux : nb_proc:%d nb_block:%d \n",size,nb_blockA);
	return 0;
      }
      nb_block = nb_blockA;
      taille_mat = taille_matA;
      buf[0] = nb_block;
      buf[1] = taille_mat;
    }
   
  //envoi de nb_block et taille_mat à tous les processus
  MPI_Bcast(buf, 2, MPI_INT, 0, MPI_COMM_WORLD);
  nb_block = buf[0];  
  taille_mat = buf[1];
  
  taille_block = taille_mat/nb_block;
  block_s = taille_block*taille_block;

  //Création de la grille de taille_block²
  int dims[2] = {nb_block,nb_block};
  int period[2] = {1,1};
  MPI_Cart_create(MPI_COMM_WORLD,2,dims,period,1,&GridComm);
  int coord[2];
  int gridrank;
  MPI_Comm_rank(GridComm,&gridrank);
  MPI_Cart_coords(GridComm,gridrank,2,coord);
 
  //allocation des blocks locaux
  // block_s : nbre d'éléments d'un bloc
  int * a = malloc(sizeof(int)*block_s);
  int * b = malloc(sizeof(int)*block_s);
  int * c = malloc(sizeof(int)*block_s);
  int z;  
  for (z=0;z<block_s;z++) c[z]=0;
  int * tempA = malloc(sizeof(int)*block_s);
  
  //allocation des matrices pour les processu autre que root
  if (myrank > 0)
    {
      A = malloc(sizeof(int)*taille_mat*taille_mat);
      B = malloc(sizeof(int)*taille_mat*taille_mat);
    }

  //Envoi des matrices entières.
  MPI_Bcast(A,taille_mat*taille_mat, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(B,taille_mat*taille_mat, MPI_INT, 0, MPI_COMM_WORLD);
  
  //Récupération des blocs intéressant:
  int k = 0;
  int r = 0;
  int j = 0;
  int start = (coord[0]*taille_mat + coord[1])*taille_block;
  while (j<taille_block*taille_block)
	{
	  for (k=0 ; k<taille_block ; ++k)
	    {
	      a[j] = A[start + k + r];
	      b[j] = B[start + k + r];
	      j++;
	    }
	  r+=taille_mat;
	}
 
  /*
  printf("a\n");
  print_block(a,taille_block);  
  printf("b\n");  
  print_block(b,taille_block);  
  */

  free(A);
  free(B);
  
  //Création InterComm en ligne
  MPI_Comm LineComm;
  MPI_Comm_split(GridComm,coord[0],gridrank,&LineComm);
  int linerank;
  MPI_Comm_rank(LineComm,&linerank);
  //printf("[%d,%d] -> %d\n",coord[0],coord[1],linerank);
  gettimeofday(&mid,NULL);
  //calcul
  int pas;

  for (pas=0 ; pas<nb_block ; pas++)
    {
      
      //Bcast du bloc A;
      if(linerank == (pas + coord[0])%nb_block)
	{
	  //Copie de a dans tempA
	  for(k=0 ; k<block_s ; ++k)
	    {tempA[k] = a[k];}
	}
      MPI_Bcast(tempA,block_s, MPI_INT, (pas + coord[0])%nb_block, LineComm);
      MPI_Barrier(GridComm);
      
      /*
      printf("tempA\n");
      print_block(tempA,taille_block);  
      */

      //calcul effectif
      
      block_mult(tempA,b,c,taille_block);
      //      printf("processus %d à l'étape %d : a=%d  b=%d tempA=%d c=%d\n",myrank,pas,a[0],b[0],tempA[0],c[0]);
      /*     
      printf("c\n");        
      print_block(c,taille_block);  
      */

      //Envoi du bloc B;
      int upperrank;
      int underrank;
      MPI_Cart_shift(GridComm,0,-1,&underrank,&upperrank);
      MPI_Sendrecv(b,block_s,MPI_INT,upperrank,tag,b,block_s,MPI_INT,underrank,tag,GridComm,&status);
      MPI_Barrier(MPI_COMM_WORLD);
      //printf("coord[%d,%d], pas = %d,b[0]=%d,b[1]=%d,b[2]=%d,b[3]=%d\n",
      //     coord[0],coord[1],pas,b[0],b[1],b[2],b[3]);
      MPI_Barrier(MPI_COMM_WORLD);
    }

  //reduction de la matrice obtenue
  int * tempC;
  int i = 0;
  if(gridrank == 0)
    {  
      tempC = malloc(sizeof(int)*taille_mat*taille_mat);
    }
  

  MPI_Gather(c,block_s,MPI_INT,tempC,block_s,MPI_INT,0,GridComm);
  MPI_Barrier(MPI_COMM_WORLD);
  
  if(gridrank == 0)
    {
      
      /*  
	  printf("tempC\n");
	  print_block(tempC,taille_block);  
      */
      
      //Recopie de la matrice dans res pour reconstituer les blocs.
      int * res = malloc(sizeof(int)*taille_mat*taille_mat);
      
      int q;
      int r;
      int col;
      int ligne;
      int ligneB;
      int colB;
      i=0;      
      while(i<taille_mat*taille_mat)
	{
	  q=i/block_s;
	  r=i%block_s;
	  //dans la matrice par blocs	  
	  col = q%nb_block;
	  ligne = q/nb_block;
	  //dans le bloc
	  colB = r%taille_block;
	  ligneB = r/taille_block;
	  
	  res[ligne*taille_block*taille_mat + col*taille_block + colB + ligneB*taille_mat] = tempC[i];
	  //  printf("res[%d] = tempC[%d];\n",ligne*taille_block*taille_mat + col*taille_block + colB + ligneB*taille_mat,i);
	  i++;
	}
      gettimeofday(&end,NULL);
      printf("temps total: %lu\n",((end.tv_sec - init.tv_sec)*1000000+(end.tv_usec-init.tv_usec)));
      printf("temps d'initialisation: %lu\n",((mid.tv_sec - init.tv_sec)*1000000+(mid.tv_usec-init.tv_usec)));
       /*printf("res"\n);
	 print_block(res,taille_block);  
  */
      //ecriture de res dans un fichier
      char * matC = argv[3];
      FILE * fd = fopen(matC,"w+");
      for (i=0 ; i<taille_mat ; ++i)
	{ 
	  for(k=0 ; k<taille_mat ; ++k)
	    {
	      fprintf(fd,"%d ",res[i*taille_mat+k]);
	    }
	  fprintf(fd,"\n");
	}
      fclose(fd);
    }

  MPI_Barrier(MPI_COMM_WORLD);

  

  free(tempA);
  free(a);
  free(b);
  free(c);
  MPI_Finalize();
  return EXIT_SUCCESS;
}
